clear;clc;
sx=[0 0.5; 0.5 0];                  % matrix of S_x for spin 1/2
sy=[0 -0.5*1i; 0.5*1i 0];           % matrix of S_y for spin 1/2
sz=[0.5 0; 0 -0.5];                 % matrix of S_z for spin 1/2
od=[1 0; 0 1];   
rho=[0 0 0 0;0 1 -1 0; 0 -1 1 0; 0 0 0 0]*1/2;
sx1=kron(sx,od);sy1=kron(sy,od);sz1=kron(sz,od);
sx2=kron(od,sx);sy2=kron(od,sy);sz2=kron(od,sz);
sx12=sx1*sx2;od2=eye(4);
ax12=trace(sx12'*rho)/norm(sx12)
aE=trace(od2'*rho)/norm(od2)